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Abstract 

We investigate the modulational instability of nonlinear Schrodinger equations with periodic variation 
of their coefficients. In particular, we focus on the case of the recently proposed, experimentally realizable 
protocol of Feshbach Resonance Management for Bose-Einstein condensates. We derive the corresponding 
linear stability equation analytically and we show that it can be reduced to a Kronig-Penney model, which 
allows the determination of the windows of instability. The results are tested numerically in the absence, 
as well as in the presence of the magnetic trapping potential. 
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1. Introduction 



The experimental realization and intense theoretical studies of Bose-Einstein condensates (BECs) £Q have led 
to an explosion of interest in the field of atomic matter waves and nonlinear excitations in them. In particular, 
one-dimensional (ID) dark [2] and bright 3 solitons have been created in experiments, in BECs with repulsive 
and attractive interactions, which correspond to the positive and negative atom-atom scattering length (SL), 
respectively. Additionally, multi-dimensional solitons 4 , Faraday waves [S], ring dark solitons and vortex 
necklaces [jj, have been theoretically predicted. The BEC solitons are of fundamental interest, not only 
conceptually, but also practically: one can envision robust solitary structures being coherently manipulated 
in matter- wave chips, similar to how light is controlled in the existing optical ones 

To generate BEC solitons in one spatial dimension (in the presence of the magnetic trap |B1 EH , or m 
the quasi-discrete setting created by an optical-lattice potential ) , as well to avoid their collapse in two 
dimensions an experimentally realizable protocol has been recently proposed, in the form of the so-called 
Feshbach Resonance Management (FRM). It is based on adjusting the effective SL (including a possibility 
to change its sign) by a resonantly tuned ac magnetic field through the Feshbach resonance |12j . 

The FRM scheme can be modelled (in the mean-field approximation) in the framework of the Gross- 
Pitaevskii (GP) equation pQ, with the coefficient in front of the nonlinear term being a periodic function of 
time. In Ref. |5], the periodic function has been taken to be a piece- wise constant one, periodically jumping 
between positive and negative values. The same model may also be realized in terms of nonlinear optics, 
where it applies to a medium composed of alternating layers with opposite signs of the Kerr nonlinearity 
|13| . FRM resembles the dispersion-management (DM) scheme, well-known in fiber optics (see, e.g., |14) 
and |15 | for review), which assumes that nonlinear fibers with opposite signs of the group- velocity dispersion 
periodically alternate, to form a system supporting robust breathing solitons. Similarly to this, very robust 
and stable matter- wave breathing solitons can be generated by means of the FRM in BECs [B], which lead 
to structures (such as, e.g., the dark solitons) even more robust than their optical counterparts 

The purpose of this paper is to study how solitary waves can arise in the FRM setting and, in particular, 
to investigate the modulational instability (MI) in the corresponding nonlinear Schrodinger (NLS) equations 
with periodically varying coefficients. Such a study is relevant, as MI is one a fundamental mechanism 
that leads to the formation of localized solitary-wave structures in a variety of settings, ranging from fluid 
dynamics (where it is usually referred to as the Benjamin-Feir instability) |17j to nonlinear optics |18| and 
plasma physics ^Hjj and, most recently, to BECs in the presence of an optical lattice fSIIS or magnetic- 
trap potentials 21 . In the present work, the study of the MI leads to a general linear stability equation 
which includes, as special cases, equations previously examined in the context of DM. We then proceed 
to analyze MI conditions in the case of the piecewise-constant FRM scheme [H], and find that the resulting 
stability equation resembles the well-known Kronig-Penney model of solid-state physics |22j . Thus, we derive 
analytical criteria for the MI. 

The paper is organized as follows. The derivation of the stability equation and its analytical treatment 
are presented in Section 2. In Section 3, we turn to numerical experiments both in the absence of the 
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magnetic trap (the case for which we have analytical results) and in the presence of it, as the magnetic 
trap is a necessary ingredient of the experiment. In Section 4 we summarize our findings and present our 
conclusions. 



2. Stability analysis 

We consider a general GP/NLS equation, with time-dependent coefficients in front of the dispersive and 
the nonlinear terms, D(t) and a(t), respectively. The equation also includes a term corresponding to the 
external magnetic (parabolic) potential trap: 

iut = ~D(t)u xx + a(t)\u\ 2 u + {l/2)n 2 x 2 u. (1) 

In the case of BECs (D = const = 1/2), u(x, t) is the macroscopic wave function, x and t are measured in 
the harmonic-oscillator units, and Vl is the strength of the magnetic trap [Q. If the FRM is applied [Hj, the 
nonlinear coefficient a(t) periodically alternates between a\ (for < t < r) and —0,2 (for r < t < T), where 
it is assumed that 0.1,02 > 0. As Eq. is spatially inhomogeneous due to the presence of the parabolic 
trap, analytical investigation of the linear stability of the uniform states is only possible for = 0. Hence, 
we will examine this case, which will be complemented by numerical results for both f2 = and Q ^ 0. 
The plane-wave solution to Eq. is 



ito = Aq exp 



i{—q 2 / D(s)ds — A / a(s)ds + qx 
Jo Jo 



Notice that the solution's amplitude can be rescaled to Aq = 1. Then, a solution including an infinitesimal 
perturbation is sought as 

u = uq [1 + ew(t) cos(fcir)] , 

where e and k are the amplitude and wavenumber of the perturbation, which leads to the following linearized 
equations for w = w r + iwi : 

w r = k 2 D{t)wi, ibi = - [k 2 D{t) + 2a(t)] w r , 
which can be transformed into a single equation, 

w r = DD- 1 ^. - k 2 D(t) [k 2 D(t) + 2a(t)] w r , (2) 

the overdot standing for d/dt. 

There are several special cases of this equation that were previously studied. In the case of DM (i.e., for 
D = D(t) and a(t) = const), the MI analysis was performed in On the other hand, in the FRM context, 
for D = 1/2 and time-periodic a(t), Eq. (j2J is a Hill equation that was considered in [H], while the more 
specialized case of a(t) = 1 + 2acos(ujt) leads to the Mathieu equation that was dealt with, in this context, 
in|S|. 

Motivated by the FRM scheme proposed in [H], we will now explore the special case of a piecewise- 
constant time-dependent SL (alternating between a\ and —02 as discussed above). In this case, we define 
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sf = fc 2 (fc 2 /4 + a,i) and s\ = k 2 (a 2 — k 2 /4), which assumes k 2 < 4a 2 . Searching for a solution in accordance 
with Bloch's theorem |22| . we find, within one period, 

w r = Ae^+^ + Be^-"-* 1 ^, for < t < r, (3) 
w r = C* e (-^+«2)t + £)e ( - iw - S3 )*, for r<t<T. (4) 

Then one should request the continuity of the solution and its derivative at t = and t = r, obtaining 4 
conditions on A, B,C, D, namely, 

e to »M + e _falT B-e" 3,r C-e _a!, ' r D = 0, 
+ si)e isiT 7l + - si)e" lslT B - (-iw + s 2 )e S2r C - (-iw - s 2 )e~ S2T D = 0, 

A + B - e (-i"+S2)T C _ e (-iu-s 2 )T D _ 0j 

+ Si)A + t(-w - - (-iw + s^e^+^C - (-tw - s 2 )e { - luJ - S2)T D = 0. 

The resulting system of homogeneous linear equations for (A, B,C, D) has non-trivial solutions only if 
its determinant vanishes. This condition leads to the following equation for the eigenfrequency (Floquet 
multiplier) ui: 

2 2 

cos(wT) = - 5 * - 52 sin(sir) sinh[s 2 (T - r)l + cos(sir) cosh[s 2 (T - r)l = F(fc). (5) 
2sis 2 



If the above condition k 2 < 4a 2 does not hold, we redefine s 2 = ^Jk 2 (k 2 /4 — a 2 ), and obtain, instead of Eq. 
©, 

2,-2 

cos(wT) = - 1 _ 2 sin(sir) sin[S 2 (T - t)1 + cos(sir) cos[S 2 (T - t)] = (6) 

By examining the function |f(fc)| or F(k) , defined in Eqs. (J5J and 0, and comparing it to 1, we 
can find whether there is a real eigenfrequency for a given perturbation wavenumber k, or it belongs to a 
"forbidden zone" , which implies the MI. 

A typical example of results produced by this analysis is shown in Fig. ^ for a\ = a 2 = 0.3 and 
r = T/4 = 1. The first several instability windows, which are determined by the shape of the curve F{k) 
for this case are given in Table I. Notice that in the table, only positive values of unstable wavenumbers are 
given, the negative ones obeying the k — > —k symmetry. In general, we have found that the number and 
widths of the instability windows increase as long as the mean value a = [a\T — a 2 (T — t)]/T of the SL gets 
large and negative, i.e., when the BEC is "very attractive" on the average. 

3. Numerical Results 

We now examine the validity of the above theoretical predictions through full numerical simulations of Eq. 
ifTJl . This also allows us to investigate the role of the magnetic trap that was not taken into regard in the 
above analysis. 

As a typical example (similar results have been obtained in other cases) , we consider modulationally stable 
and unstable cases for the example presented above (i.e., a\ — a 2 = 0.3 and r = T/4 = 1.). In particular, 
we examine the case of k — 0.5 which, according to Fig. n an d Table I, should be unstable, and the case 
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of k = 1, that is expected to be stable. It can be clearly observed that the corresponding configurations 
evolve in accordance with the theoretical prediction, indeed leading to the MI and, as a result, generation 
of a lattice of solitary waves, for k = 0.5. Instead, the modulational perturbations remain small for k = 1. 
The agreement between the analytical predictions and direct simulations has been found to be generic. 

We now proceed to the case when a weak magnetic trap is included, taking (for example) £1 2 = 0.00025. 
Then we observe that both the configurations which were linearly stable and unstable in the absence of the 
trap eventually develop the MI, which is in consonance also with the findings of |21| . Here, we used the initial 
condition u = utf[1 + 0.02 cos(fc#)], where utf == \/max(0, /i — V(x)) is the wave function in the Thomas- 
Fermi approximation £Q. The reason why the MI develops in both cases is that the magnetic trap mixes 
stable and unstable wavenumbers in the Fourier space. As a result, even though the initial configuration 
contained only stable wavenumbers, unstable ones are eventually generated, giving rise the MI. However, as 
can be observed in Fig. the MI occurs faster in the case that would be modulationally unstable in the 
absence of the trap, as the MI growth rate is larger in this case than in the case when the uniform state was 
stable in the absence of the trap. 

We also examined what happens with the increase of the strength of the magnetic trap, which makes it 
"tighter" . In this case, the trapping potential is much more efficient in mixing the wavenumbers, therefore 
there remains little difference between the evolution of configurations that were modulationally unstable and 
stable without the trap. An example is shown in Fig. 0]for Q 2 = 0.0015 (6 times as large as in the case 
shown in Fig. EJ). It can be seen that the development of the MI for k = 0.5 and k = 1 is practically identical. 

4. Conclusions 

In this paper, we have studied the modulational stability (MI) in the GP/NLS equations with periodically 
varying coefficients. We highlighted how the results may be applied to problems stemming from both 
nonlinear optics and Bose-Einstein condensation. Findings reported in earlier works in the subject have 
been generalized. The equation obtained from the analysis, Eq. J2J), is of interest in its own right. Motivated 
by the recently introduced concept of a periodically modulated (via a Feshbach resonance) scattering length 
as a means of generating robust breathing matter- wave solitons [B] , we have focused attention on the case 
of alternating positive and negative scattering length. 

For the case of the FRM, we have demonstrated that the above-mentioned equation governing the MI 
development is an equation of the Hill type. More specifically, it is an analog of the Kronig-Pcnney model 
known in solid state physics. Then, developing analysis similar to the eigenvalue calculation in the Kronig- 
Penney model, we have obtained analytical conditions for the MI of such FRM schemes. The analytical 
criteria have been shown to be in agreement with full numerical simulations. Although the analysis was 
carried out in the absence of the magnetic trap, the role of the trap was investigated numerically. For weak 
trap potentials, it was found that the MI will always set in (due to the mixing of wavenumbers imposed by 
the spatial potential); however, its growth rate bears "memory" of the presence/ absence of the MI in the 
corresponding situation without the trap. On the other hand, if the trap is strong, the memory is effectively 
wiped out. 
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Window 


k 


1 


[0,0.7489) 


2 


(1.2502,1.4957) 


3 


(1.8283,1.9115) 


4 


(2.2396,2.2562) 


5 


(2.5638,2.5784) 



Table I: The first instability bands of the perturbation wavenumber k for the case a\ = a 2 = 0.3 and 
r =T/4= 1. 




k 



Figure 1: The plot of the function F(k) defined by Eqs. © and © for ai = a 2 = 0.3 and t = T/4 = 1. 
When the function satisfies \F (k)\ < 1, the perturbation of wavenumber k is modulationally stable, while 
for |-F(fc)| > 1, it is modulationally unstable. 
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Figure 2: For the case a\ = a 2 = 0.3 and r = T/4 = 1, fi = 0, the evolution of modulationally unstable 
(k = 0.5, left panel) and stable (k = 1, right panel) perturbations of the form 0.02 cos(kx) added to u = 1, is 
examined. The top and bottom panels show, respectively, the configuration at t = 40 and the time evolution 
of the amplitude of the spatial profile. 





20 

t 




Figure 3: The same as in the previous figure, but now with f2 2 = 0.00025, the initial configuration being 
multiplied by the Thomas- Fermi wave function ^/max(0, [i — V{x)) with /i = 1 . Notice that the modulational 
instability develops both for k = 0.5 (left) and for k = 1 (right). 
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Figure 4: The same as in the previous figure, but with fi 2 = 0.0015. The only remaining visible difference 
between the cases of k = 0.5 (left) and k = 1 (right) is a slightly larger growth rate in the former case. 
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